Make sure you have installed the latest version of the connectome workbench, which you can download here: https://www.humanconnectome.org/software/connectome-workbench
#### load required libaries ####
require(neurobase)
require(ciftiTools)
require(plyr)
require(dplyr)
#### set path to workbench (example for a MacOS environment) ####
ciftiTools.setOption('wb_path', '/Applications/workbench/bin_macosx64')
## Using this Workbench path: '/Applications/workbench/bin_macosx64/wb_command'.
For this example, we will extract ROI based data from a tau-PET image using the 17 network version of the Schaefer 200 ROI atlas. The atlas including the reference (Schaefer et al., Cerebral Cortex, 2018) can be found here: https://github.com/ThomasYeoLab/CBIG/tree/master/stable_projects/brain_parcellation/Schaefer2018_LocalGlobal
ste the path to your downloaded tutorial
sdir = "/Users/isd_neuroimaging/Library/Mobile Documents/com~apple~CloudDocs/Github/Tutorials/Connectome_workbench/Map_nifti_to_ROIs/"
# load the atlas
atlas_readin <- neurobase::readnii(
file.path(sdir, "/Example_data/rSchaefer2018_200Parcels_17Networks_order.nii.gz")
)
neurobase::ortho2(atlas_readin)
# load the atlas
pet_readin <- neurobase::readnii(
file.path(sdir, "/Example_data/smoothed_8mmFWHMwsub-0555_ses-2020-02-06_pet-AV1451_flirtT1reg_INFCER_SUVR.nii.gz")
)
## Warning in datatyper(nim, warn = warn): Need to change bitpix and datatype to
## FLOAT64 due to NAs
neurobase::ortho2(pet_readin)
# check overlay
neurobase::ortho2(pet_readin, atlas_readin)
# check if dimensions match
dim(pet_readin) == dim(atlas_readin)
## [1] TRUE TRUE TRUE
pet_atlas_based =
data.frame(pet = as.vector(pet_readin@.Data),
atlas_idx = as.vector(atlas_readin@.Data)) %>%
filter(atlas_idx !=0) %>%
group_by(atlas_idx) %>%
summarize(mean = mean(pet))
head(pet_atlas_based)
## # A tibble: 6 × 2
## atlas_idx mean
## <int> <dbl>
## 1 1 1.11
## 2 2 0.742
## 3 3 0.860
## 4 4 0.904
## 5 5 1.09
## 6 6 1.16
# show distribution of SUVR data
hist(pet_atlas_based$mean)
# read in atlas
atlas_dscalar <- read_xifti(file.path(sdir, "/Example_data/Schaefer2018_200Parcels_17Networks_order.dscalar.nii"))
atlas_dscalar_vec <- c(as.matrix(atlas_dscalar))
# determine values to map to surface
vector_to_map = pet_atlas_based$mean
# map ROI-based PET data to cifti
xii_to_map <- c(NA, vector_to_map)[atlas_dscalar_vec + 1]
xii1 <- select_xifti(atlas_dscalar, 1)
xii_to_map <- newdata_xifti(xii1, xii_to_map)
# view cifti directly in R
view_cifti_surface(xii_to_map, colors = "magma", zlim = range(vector_to_map),
borders = T)